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' Abstract. Many stochastic time series can be described by a Langevin equation composed of a deterministic and 

a stochastic dynamical part. Such a stochastic process can be reconstructed by means of a recently introduced 
nonparametric method, thus increasing the predictability, i.e. knowledge of the macroscopic drift and the microscopic 
diffusion functions. If the measurement of a stochastic process is affected by additional strong measurement noise, the 
reconstruction process cannot be applied. Here, we present a method for the reconstruction of stochastic processes in 



> 

C/3 ■ the presence of strong measurement noise, based on a suitably parametrized ansatz. At the core of the process is the 

, ^ ' minimization of the functional distance between terms containing the conditional moments taken from measurement 

C/3 , data, and the corresponding ansatz functions. It is shown that a minimization of the distance by means of a simulated 

annealing procedure yields better results than a previously used Levenberg-Marquardt algorithm, which permits a rapid 
' and reliable reconstruction of the stochastic process. 



> 

^ ', 1. Introduction 

Physical systems often can be resolved on two different time scales, observing slowly varying 
^ , macroscopic motion while much faster microscopic interactions are perceived as an effective noisy 

^ I driving. An adequate description of these systems can be obtained by a Langevin equation, which yields 

a deterministic drift term and a stochastic diffusion term. With respect to this description, reconstructing 
^ I the drift and diffusion terms of an unknown process in time or scale from a data set enables one to 

enhance the predictability of that process. Recently, a nonparametric method of reconstruction based 
solely in the measured time series has been suggested [1]. The method has been applied successfully to 
^ I data sets from fields such diverse as turbulent fluid dynamics [2i, human movement [31> financial data 

^ ' S, climate indices 151161, and to electroencephalographic recordings from epilepsy patients ||7]|8]. For 

- ■ ■' the case of finite time steps ||9l[T0l, which can occur when the data is not available at sufficiently high 

sampling rates, recent improvements have been suggested lllTl[T2l[T3l . 

However, any experimental setup will record additional measurement noise when recording a data 
series, which cannot be easily distinguished from the intrinsic dynamical noise. The presence of 
measurement noise can destroy the Markovian properties of the data and lead to incorrect estimates 
for the drift and diffusion temis llT8]l20l. 

It has been suggested to separate the measurement noise from the dynamics of the measured variable 
using different predictor models or schemes for noise reduction |[T4l [T5l . Recently, another approach 
has been suggested, which allows reconstructing the properties of stochastic time series affected by 
exponentially correlated Gaussian noise, using algebraic properties of the noniial distribution ||T61 . Here, 
we revisit an approach presented in Refs. lITTl ITSl . which minimizes the functional distance between 



characteristic functions extracted from the measurement data and corresponding terms calculated from 
a suitably parametrized ansatz for the drift and diffusion functions. We show that using a Simulated 
Annealing (SA) |[T9l method for the optimization yields results that are superior to the ones obtained 
previously by means of the Levenberg-Marquardt (LM) method. 

We start in Sec. |2]by describing the general framework for extracting the evolution equation of a 
signal spoiled by strong measurement noise, together with the amplitude of said measurement noise. In 
Sec.[3]both LM and SA methods are addressed and compared. Section |4] concludes the paper. 



2. Langevin approach for data sets with measurement noise 

We consider a one-dimensional Langevin process x{t) defined as 

dx 



— = Di{x) + .jD2{x)Tt, (1) 

where Tt represents a Gaussian (5-correlated white noise (Tt) = and (^t^t') = 2(5(t — t'). The 
Kramers-Moyal (KM) functions Di{x) and D2{x) are defined as 

Dn{x) = - lim-Mn{x,T) (2) 

for n = 1,2, where Di describes the deteiTninistic drift and D2 corresponds to a diffusion process. An 
estimate of the n-th order conditional moments of the data M„(x, r) can be extracted directly from the 
measured time series as 

Mn{x^,T) = {{x{t + t) - x{t)r)U(^t)=x, , (3) 

the hat indicating that they are calculated from the measured data x{t) directly, whereas the brackets <> 
denote averaging over suitable parts of the time series - or ensemble averages, if available - as described 
in Refs. ||Il|4l|6l|T3. Thus, a nonparametric reconstruction of the stochastic process can be achieved. 

In the following, we will consider a time series generated by integrating Eq. ([T|l with drift and 
diffusion coefficient assumed to be linear and quadratic forms respectively, 

Di{x) = dio + diix (4a) 



D2{x) = d2Q + d2lX + d22X^ ■ (4b) 

Though we concentrate on the particular expressions for Di and D2 given above, it should be stressed 
that they comprehend a large collection of different processes, such as Omstein-Uhlenbeck processes 
ifTTl . The parameter diQ can always be eliminated through a simple transformation x — > = 
X + dio/dii. For the numerical examples in the remainder of the text, we will use the parameters 
diQ = = -1,^20 = l,c?2i = -1,^22 = 1- 

We now consider the case that the signal of x{t) is affected by a Gaussian 5-coiTelated measurement 
white noise, which leads to the series of observations 



y{t) = x{t) + aC{t) (5) 

where a denotes the amplitude of the measurement noise. As Figs. [T^ and [IJ5 show, although the trend of 
x{t) is still visible in the presence of measurement noise, it is intuitively difficult to distinguish between 
the intrinsic dynamic noise and the measurement noise. Likewise, Fig. [1^ illustrates that for increasing 
a one obtains broader probability density functions P{y): the standard deviation of these distribution 
increases with the measurement noise, whereas the mean value remains constant. 

In the presence of measurement noise, the conditional moments yield non-zero finite values at 
r = and thus the limit lim -Mnix^a / 0, r) does not exist lITSl . Likewise, the requirements for 

Markovianity of the measured time series y{t) may not be satisfied ll20l . 




Figure 1. Langevin time series with different measurement noise strengths. Here we show a time series 
(a) without measurement noise, and (b) with strong measurement noise. In (c) the probability density 
function -P(y) of the series with measurement noise (see Eq. dD) is shown. In all cases, the assumed time 



series x{t) without measurement noise uses the coefficients Di{ 
when integrating Eq. ([T). 



1 — X and D2{x) = 1 — x + x 



However, the two conditional moments H] |4l |6l [T71 [181 Mn{yi,T) (n = 1, 2) can be derived taking 
y{t) instead of x{t) in Eq. Q. Figure |2] shows plots of both conditional moments as a function of r for 
different strengths of the measurement noise. The linear dependence of Mi and M2 on r makes it still 
possible to obtain an estimate for the 'spoiled' KM-coefficients: 



Dniy) 



Mniy,T2) - Mniy,n] 

n\{T2 - Ti) 



(6) 



Within this description, an estimate for the amplitude of the measurement noise yields (J ^ \j ^"^^^^'^^^ 
(see Ref. ifTTl 1211 '). where /x is the average value of y{t) data points in the time series. As shown in 
Fig.|2t, such an estimate is not valid for sufficiently strong measurement noise, i. e. cj > 0.5. 

The coefficients Di and D2 cannot be correctly estimated, though. Figure |3]shows how the estimated 
parameters dij deviate from the 'true' values dij of the process given by Eq. as measurement noise 
increases. Notice that for a = (left vertical axis in each plot of Fig. O the estimated parameter values 
are approximately correct. 

To correctly derive the drift and diffusion coefficients Di {x) and D2 (x) when a is strong, we consider 
the measured conditional moments Mi{yi,T) and M2(yi,r). Since these conditional moments depend 




Figure 2. Conditional moments (a) Mi (y^ , r) and (b) M2 {ui , r) as a function of r, for bin Ui = and 
different measurement noise strengths a. The same x{t) as in Fig. [T]was used, (c) First estimate of the 
measurement noise (see text). 




Figure 3. Noise dependence of functions Di{y) and 1)2 (y) (see text and Eq. Q). The underlying 
Langevin time series x{t) without noise is the same as in Fig.[T] 



in a non-trivial way on both time r and amplitude yi, we approximate them up to first order on r: 

Mi(?/„r) = {y{t + T)-y{t))\y^t)=y^=Tmi{yi)+^i{y.i)+0{T^), (7a) 

M2{y^,T) = {{y{t + T)-y{t)f)\y(t)=y,=Trn2{y^)+^2{yi) + (y^ + 0{t''), (7b) 

where y{t) is taken in the range yi it Ay/2 for each bin i, and Ay depends on the binning considered. 

While the functions rhj and 7j {i = 1, 2) are obtained explicitly for each bin value yi, the decisive 
aspect of this approach is the fact that the ansatz functions nii and 7^ can be calculated from the drift and 



diffusion coefficients Di , D2 and the measurement noise distribution as follows lITSll : 

r+co 

li{y) = / [x -y)fa{x\y)dx (8a) 



72 (y) = / {x-yYUx\y)dx (8b) 

/+00 
Di{x)f„{x\y)dx (8c) 
-00 

/4-00 
[{x - y)Di{x) + D2{x)]f„{x\y)dx, 
-00 

(8d) 

where fa{x\y) is the probability for the system to adopt the value x if a measured value y is observed. 

The problem we then solve is to pai^ametrize Di,D2 and the noise strength a and to determine the set 
of parameters that minimize the functional distance F between the measured functions rrij, 7^ and their 
counterparts ?nj, 7^ defined as 



^' ^{ii -ii{yi)f {12 - i2{yi) - (y'^f 
.2 + ^2 ' + 



_1 V- r (7i - 

M a 



71 ' 72 * 

(9) 



(rhi - mi{yi)f ^ (rfi2 - m2{yi)f 



where the summation extends over all M bins, and ax^^{yi) is the error associated to function 71 at the 
value yi (and similarly for ax^^, a^i and (1^2, all of them taken directly from the data only). 

After computing the functions 71, 72, we start from the initially estimated set of values for the 
parameters and iteratively improve the solution by seeking lower values of F, until convergence is 
attained. In the following section we use two different methods for minimizing the functional F, namely 
the LM method f\M and SA method. We will show that for our purposes SA is significantly better than 
LM. 



3. Comparing two optimization methods: Levenberg-Marquardt and Simulated Annealing 

In this section, we consider the parameters cr, dn, 1^20, c?2i and ^22 which we denoted as pk with 
k = 1, . . . , 5, respectively. For the Levenberg-Marquardt procedure one computes the derivatives of 
F with respect to the parameters p = {pk}, and uses them to define the coefficients /3fc = ^■§^ and 

aki = Ej=i -^sj^ ^^'dpf^ ^^'dpi^^ ^ where fj{yi, p) is the j-th summand in the bracketed inner 

sum defining F, Eq. The descent 5pi is computed from X^/^j^ ^'ki^Pi — Pk^ which uses a damped 
version of the Gaussian matrix: 

/ { a'ki = {l-X)akuk = l 
^ki = \ , , ^ , (10) 

The parameter A is updated during the optimization. The procedure stops after sufficient convergence 
has been achieved. For details, see Ref. ll22l . We generated several data sets from Eq. ([1} with different 
'input' measurement noise amplitudes < o"/ < 1.2, in order to compare the reconstructed parameters 
to the actual values. 

The optimization results are shown in Fig. |4] triangles indicating the first estimate for the parameter 
values (see Sec.|2l), solid lines indicating the true values used to generate the data, and squares indicating 
the value after LM-optimization. 
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Figure 4. Comparison of the optimized parameters values (bullet) with the first estimate and the true 
values for different input measurement noise strengths aj: (a) 2cr^, (b) dio, (c) du, (d) ^20, (e) ^21, 
(f) ^22- The measurement noise is correctly extracted, as well as the parameters representing the drift 
coefficient Di{x) and the diffusion coefficient D2{x). SA results are significantly better than LM results 
(see text). Error bars represent the maximum and minimum values encountered in 10 realizations. 



While the parameters are well estimated, particularly the magnitude of measurement noise and the 
drift parameters dio and du, it was found that LM sometimes does not converge, or converges to a local 
minimum, the latter fact becoming visible when comparing the results after starting from different initial 
values. Both these shortcomings required that the LM optimization results were drawn as the best results 
of a series of attempts with different starting values. The intent to overcome these shortcomings led us 
to attempt the SA method (squares in Fig.|4]l. 

Simulated annealing is a probabilistic global optimization method. Starting from initial guesses for 
the parameter vector it proceeds with a step in a random direction. The step is accepted immediately if 
the energy function at the new position .Fnew (see Eq. Q) is lower than the previous energy Fold, or it 
is accepted with a probability given by a Boltzmann factor exp(— AF/fcT), where AF = Fncw — -^old- 
Otherwise the step is refused and a new parameter step is chosen. This procedure corresponds to the 
motion of thermal atoms in an attractive potential. Gradual cooling is then achieved by reducing the 
temperature parameter, thus annealing the test particle in the minimum of the energy landscape. 

Figure H] shows that both optimization methods determine the value of input measurement noise fj/ 
correctly in all cases, which we find remarkable given the comparable magnitude of measurement noise 
and the time series x{t) for cr « 1. Furthermore, it can be seen that both LM and SA estimate the 
parameters dik for the drift and d2i for the diffusion coefficient correctly; as expected lITSl . the drift 
coefficients are determined with higher precision. In general, the SA routine provides more accurate 
results than LM. 

In the LM case, for stronger measurement noise amplitudes, namely for a > 1.2, the algorithm is 
sometimes stuck in a local minimum of the function F, leading to unreliable coefficients djjf The same 
set of 10 time series with 10^ data points were used for each method, but in the LM case, not all cases 
showed local minima sufficiently close to the true parameter values. Therefore, in our comparison, we 
overestimate the best performance of the LM method, by considering only the best five cases out of the 
set of 10 time series (squares in Fig. |4]). 
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Figure 5. (a) Plot of the relative energy in comparison to the true value, Eq. dTTI) . as a function of the 
Euclidean distance dr = \\pi — Pi.o|j2 of the parameters from the true parameter values {pifi} found 
with SA method (circles) and with the LM method (squares) for a = 1.2. (b) Close-up of (a) near the 
minimum of F. For better visibility, not all SA-points are shown in (a) and (b). (c) Execution time (in 
seconds) as a function of a for LM (squares) and SA (circles). 



Figure |5^ shows a plot of the energy, i.e. the distance 

dF = Fip,) - Fo{p^,o) (11) 

where Fo(Pi,o) is the value found when taking the true parameter values, as a function of the Euclidean 
distance dr = \ \pi — Pifi | b of the parameters from the true parameters values. Apparently, the energy is 
strongly asymmetrical in at least one of the parameters, which is in agreement with Ref. lITSl . Also, the 
existence of local minima, which apparently attract the LM algorithm, is confirmed. Figure |5j) shows a 
close-up near the lowest minimum to emphasize the higher accuracy of the final result for SA (circles) 
compared with LM (squares). In this example, it can be seen that the LM algorithm does not converge to 
the global minimum but deviates shortly before. 

Another argument in favor of the SA algorithm is that the execution time approximately 100 s lower 
than for the LM method, as shown in Fig. |5]:. Additionally, whereas the LM algorithm needs to execute 
on 10 reahzations of the time series in order to produce a reasonably accurate result, the SA algorithm 
was found to converge to the global minimum in almost all instances, implying that fewer realizations 
are needed for averaging. 

4. Discussion and Conclusions 

We describe the improvement of a nonparametric procedure to extract measurement noise in empirical 
stochastic series with strong measurement noise by applying a simulated annealing approach for the 
optimization step. Simulated annealing accurately extracts the strength of measurement noise and the 
values of the parameters representing drift and diffusion. These parameters fully describe the evolution 
equation for the measured quantity. SA produces more reliable and accurate results than a previously 
used Levenberg-Marquardt algorithm, with the additional benefit of being significantly faster. 

Nonparametric reconstruction of the Langevin Eq. ([T) from measured stationary data sets merely 
requires that the process exhibits Markovian properties and fulfills the Pawula theorem f6l, whereas 
the measurement noise needs to be uncorrelated. The Pawula constraint can be relaxed extending the 
analysis to Levy noise |[23ll24]| . It should be pointed out that the method presented here relies solely on 



the Markovian properties of the underlying, undisturbed process x{t), and does not require the measured 
process y{t) to be Markovian. 

An extension of this work to multidimensional Langevin time series will be essential for assessing the 
complexity behind EEG time series or earthquake data. It is likely that we may combine this denoising 
approach with the method of eigendirections presented in ll25l . 
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